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Three computer programs were written to find the eigenvalues 
and eigenfunctions of acoustic normal modes in the ocean. The 
programs used two different methods: an iterative finite difference 
scheme, and a method based upon the WKB approximation of quantum 
mechanics. The methods assume a flat fluid bottom and are designed 
for any arbitrary sound speed profile. While the results of both the 
finite difference and the WKB methods agreed, the WKB method 


proved faster. 
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TABLE OF SYMBOLS 


The following symbols are used within the text of this thesis, 


but do not necessarily apply to the computer programs. The figures 


in parentheses refer to the equation which defines the symbol, or in 


which the symbol is first mentioned. - 

a coefficient for particular solution of Zy (Stk 

Ai Airy phase solution. 

b coefficient for particular solution of Z) (57). 

Bi Airy phase solution. 

Cc sound velocity. 

oe bottom sound velocity. 

Sain minimum sound velocity in the sound velocity profile. 

C COcimiclient tor particular solution of Z5 (24). 

Cy group speed (40). 

os phase speed (39). 

D coefficient for particular solution of Z5 (24). 

e 2.718281828459045. 

E energy level, or horizontal wave function (17). 

F(z) square of the vertical wavenumber (92). 

h depth grid spacing. 

H depth from upper (lower) turning point to surface (bottom) 
(OZ): 

lee A) Hankel functions of the first and second kind. 





wave number (14). 

maximum possible real wave number (16). 
horizontal wave number (14). 

eigenvalue (horizontal wave number). 
vertical wave number (14). 

bottom imaginary vertical wave number (24). 


imaginary vertical wave number at water side of bottom 
(L02.).s 


imaginary vertical wave number (49). 
integral of imaginary wave number across a barrier (86). 
mode number. 

maximum number of discrete modes available (33). 
acoustic pressure. 

horizontal range. 

horizontal displacement potential (10). 

reflection coefficient for a barrier (85). 

an integral function of depth (34). 

vertical wave number integral for Z, (65). 

vertical wave number integral for Zo (64). 

time. 

velocity function (16). 


bottom value of the velocity function. 





2b 


oor 1, 


ae 


5a 


depth. 

bottom depth. 

lower turning point depth (64). 

upper turning point depth (66). 

top depth of a barrier (86). 

bottom depth of a barrier (86). 

depth dependent vertical displacement potential (13). 
bottom vertical displacement potential (22). 


vertical displacement potential at water side of bottom 
boundary. 


eigenfunction for a particular mode, m. 
WKB solution in the ''real'' region (57). 
WKB solution in the "imaginary" region (59). 
mode attenuation coefficient (45). 

bottom volume attenuation coefficient (45). 
difference function (80). 

an undefined function of depth (34). 

phase at lower turning point (63). 

phase at upper turning point (67). 

phase at top of barrier (115). 

phase at bottom of barrier (115). 

complex wave number in the etter (43). 
complex horizontal wave number (44). 
mode normalization integral (41). 


an undefined function of depth (35). 





P, 


Co 


Q 


c OH Ss 


3. 141592653589793 
dens tive 

bottom density. 
ocean water density. 
velocity weighted mode normalization integral (42). 
angle made by ray with the horizontal (grazing). 
displacement potential (2). 


angular frequency. 


10 









I. INT RODUCIIeN 


Since the Second World War considerable interest has developed, 
in both naval and scientific communities, in the prediction of sound 
propagation within the ocean, For the scientist, propagation informa- 
tion can be a vital part in investigating the acoustic, physical and 
sometimes chemical properties of the ocean. For the naval commu- 
nity, propagation is a vital element in the prediction of acoustic sen- 
sor performance. The prediction of such performance is not limited 
to the design and development of sensor systems, but is increasingly 
important in the operational employment of such sensors and selec- 
tion of the most fruitful tactics. 

Accordingly, there has been considerable development within the 
naval community of numerical acoustic propagation models. These 
models have, for the most part, been based upon the theory of ray 
acoustics. The techniques employed have developed considerable 
sophistication in order to deal with ray limitations with caustics, 
frequency effects, and interference. However, for prediction of prop- 
agation over great ranges and at low frequencies, ray techniques 
have two serious limitations: 

(1) The long ranges involved require long computation times. 
Each “a must be ''traced'' over many successive short time steps to 


simulate travel of the wavefront. Usually more than one hundred such 


rays are traced. Thus, as ranges increase, the computation time 
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increases proportionately. Such long computation tirnes make the 
routine operational use of ray trace programs covering paths of 
thousands of kilometers too expensive to be practical. 

(2) Ray trace programs largely ignore such frequency dependent 
effects as diffraction and interference. The assumptions of ray a- 
coustics require that the wavelength of the acoustic signal be short 
in comparison to the depth and velocity gradient scale. As lower 
frequency signals are considered, the wavelength becomes an appre- 
ciable fraction of the depth scale, and so violates a basic validity 
assumption for ray acoustics. 

Because of these limitations there has been increased interest 
in normal mode theory as a basis for long range propagation models. 
The theory of normal mode propagation is not new and has had a num- 
ber of shallow water applications [Officer 1958, Bucker and Morris 
1965]. A. O. Williams (1970) has outlined a method whereby the 
normal mode technique could be applied to deep ocean acoustics. As 
a part of any such technique both the characteristic horizontal wave- 
number (eigenvalue) Pr the mode shape (eigenfunction) must be 
solved for ite mode. Todothis, three basic methods have been 
used: 

(1) The wave equation has been solved in closed form [Tolstoy 
and Clay 1966, Williams and Horne 1967]. in order to do this, the 
sound velocity profile must be fitted to an assumed analytic function. 


Such functions have taken the form of Epstein profiles [Bucker and 


eZ 





Morris 1967], profiles with a constant gradient of the reciprocal of 
the sound velocity or sound velocity squared [Williams and Horne 
1967, Tolstoy and Clay 1966]. The requirement that the sound veloc- 
ity profile be fitted to an arbitrary function places severe limitations 
on the flexibility of any model. 

(2) The wave equation has been vertically integrated while meet- 
ing appropriate boundary conditions. This method has the advantage 
of being adaptable to any arbitrary sound velocity profile. Kanabis 
(1972) and Newman and Ingenito (1972) have developed two such shal- 
low water normal mode models. The first program presented with 
this thesis, NORMOI, is based upon these two shallow water models. 

(3) The Wentzel-Kramers-Brillouin (WKB) approximation of 
quantum mechanics provides an analytic solution to the wave equation. 
This solution is singular at depths which correspond to ray vertex 
depths. However, the WKB approximation may lend itself to solu- 
tion of the eigenvalue, after which the eigenfunction may be solved 
by an integration similar to the previous method. The last two pro- 
grams presented with ee thesis, NORMO4 and NORMOS5, use this 
method. | 

The long term goal of this research is to develop a deep ocean 
normal mode acoustic propagation model of sufficient computational 
speed to be considered for operational use. In order to achieve this 
goal, a rapid method of solving for the eigenfunction must first be 
developed. Based upon this need the immediate objectives of this 


thesis are: 


lie 


(1) To compare the results of the programs employing methods 
(2) and (3) above with analytic solutions, the published results of 


other programs, and each other. 


(2) Determine whether the WKB method offers promise of im- 


frovement in terms of time versus accuracy. 


(3) To outline future improvements to either method based upon 


tie preliminary results. 
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li. TeeORY 


hee WAVE, THEORY 

In this section we apply a mathematical separation of variables 
to the wave equation; then we impose boundary conditions on the re- 
sulting separated vertical equation. The results of this mathematical 
procedure offer us a differential equation and boundary conditions 
capable of numerical solution. We will close this section with a dis- 
cussion of some general properties of such a solution. 

The following discussion is based upon the treatments given by 
A. O. Williams (1970), and I. Tolstoy and C. S. Clay (1966). Fora 
more complete descriptive treatment the reader is referred to the 
discussion given by Williams. 

1. The Wave Equation 

The simple scalar wave equation for small amplitude waves 


ee 2: 


a 
: V*o = Got 





nN 


(1) 


where 6 is the displacement potential. Here 4d’, the displacement 


of a fluid particle from its rest position, is represented by 


a= VO, @) 
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and the acoustic pressure, p, is given by 


ae eee 3) 


Thus, if we can express © as a function of space, we can then find 
the resulting acoustic pressure and sound propagation loss due to 
spreading and refraction. 

If we use a cylindrical coordinate system and assume ® is 


a function of range, depth, and time, we can rewrite (1) as 


5 a § xd 1 @& 
S-(r3e) +28 -4 F8 -o. (4) 





Let us further assume that © is separable in terms of r, z, andt, 


such that 


D(r,2.6) = Ror) Zz) a (5) 


where w is the source angular frequency, Z(z) is the vertical dis- 
placement potential function, and R(r) is the radial displacement 


potential function. We immediately notice that 


ad wg. 
| ia 7 


Using equations (5) and (6) we rewrite equation (4) as 


oe nak | da Oe) ore ; 7 
an ar ve Ce 2 ” 
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This expression has a range dependent term, a depth depen- 
dent term and aterm with c, the sound velocity. If we could assume 
that c is independent of either depth or range, the expression is sep- 
arable into depth and range differential equations. The obvious choice 
is to assume that c varies only with depth. This assumption, termed 
horizontal stratification, although not true over great distance, is 
commonly made in oceanography and ocean acoustics. Actually, for 
our purposes it will be sufficient if c is horizontally stratified in the 
vicinity of our solution, and further if any horizontal variation of c is 
small with respect to ‘Ne vertical variation. Making such an assump- 
tion, we continue by separating (7) into range and depth differential 


equations 


44a oR\ _ 32 
PR 22(r on) =-kp ; (3) 
A De fe oe = 6 

oA +e ke. (9) 


Solutions to equation (8) include the Hankel functions of first 


and second kind representing cylindrical waves, 


Remy = Fle” (her), (10) 


mor distances such that k,r is much greater than one, this function 


can be approximated by its asymptotic form 


pons 2 ti (hpi 2) 
Rn) ee e fe (11) 


ee 





We now have an assumed time dependence and a solution to the range 


dependent function of 6. Thus, at long range from the source we 


_ 2. -i (Ayr - T-ot) 
— (z 4 
P= % jae e ee 


Here k, , an eigenvalue which appeared in the separation process as 


have for > 





a mathematical constant, will be seen to be a horizontal wave number. 
We will consider only positive values of k,, which represent outgoing 
wavefronts. It remains for us to find a method of solving for Z(z) in 
terms of the separation parameter, k, 


We can rewrite (9) as 


22 -(S-#E)zZ-0. (13) 


> =* ce 


This equation is the separated, space form of the wave equation, or 
Helmholtz equation, for the vertical wave component. The expression 
in parentheses represents a vertical wave number, ko and is related 


to the (true) wave number, k, and horizontal wave number, kL» by 


he os ee (14) 


Equation (13) with a pair of associated homogeneous boundary condi- 
tions forms a Sturm-Liouville problem, which is solved in terms ofa 
set of eigenvalues and corresponding eigenfunctions. It now remains 
for us to define and employ boundary conditions associated with the 
ocean vertical sound profile. Before so doing, however, we will di- 


gress a moment to introduce a notational convenience. 
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Since the value of c does not vary greatly (perhaps five per 
cent) throughout the depth profile, the square of the vertical wave 
number represents a small difference between two nearly equal num- 


bers 


he SRE (15) 


In order to conserve accuracy and facilitate computational speed, a 
notational convention is borrowed from quantum mechanics. The 
square of the (true) wave number, ee is subtracted from an arbitrary 
constant of the same order of magnitude (here the maximum possible 


wave number) to form a "potential function." 


uy Ode 2 ode 
Ye) <a -KO- Ss (16) 


The square of the horizontal wave number (our separation constant) 
is also subtracted from the maximum wave number to form a quantity 
corresponding to the "energy level'' of quantum mechanics, and our 


new eigenvalue 


er. © he - (17) 


max 


With this notation, the Helmholtz equation (13) becomes 


FZ -[e-Va|Z=0 , een 


oz 


a form similar to Schroedinger's equation. By this notation we hope 
to facilitate the adaptation of the results of quantum mechanics to the 


problem at hand. 
19 





2. Boundary Conditions 

We continue by developing appropriate vertical boundary 
conditions. 

Figure 1 represents a typical deep ocean sound profile. 
Mreure Z is the same profile represented in terms of V(z). Note that 
at some intermediate depth the value of sound velocity reaches a min- 
imum at the sound channel axis. Note also that the function V(z) 
forms a ''potential well'' with depth, representing the deep sound 
channel. Our domain of interest is bounded by the sea surface anda 
flat bottom. 

The air-water boundary at the surface approximates a free 
surface at which the stresses vanish. Ina fluid this boundary corre- 


sponds to 
Z(o) = 0. (19) 


If the limit of the second derivative of the vertical displacement func- 
tion is finite at the surface, the stipulation that V(z) discontinuously 
approaches infinity at the surface has the effect of making Z(z) vanish 
at the ne. Thus, in terms of our quantum mechanics analogy, 
the free surface boundary condition corresponds to a perfectly rigid 
wall of the potential well. 

The ocean bottom provides a boundary more difficult to 


specify. The ocean bottom has a complex, often layered, horizontal- 


ly varying structure. In addition, little is known of the specific 
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acoustic properties of much of the bottom required tomodel fully the 
bottom boundary conditions. In the face of such complexities we will 
approximate the bottom with a relatively simple model. 

We will assume that the bottom is represented by a homo- 
geneous fluid of infinite depth, and constant velocity. Across a fluid- 
fluid interface we require that both pressure and vertical particle 
motion be continuous. From (3) and (6), continuity of pressure re- 


quires 
e, Zo = erhe (20) 


at the bottom. From (2) we find that continuity of vertical particle 


motion or displacement at the bottom requires 


ev 
NS 
) 

rw 





| 


Ze 
oz ey 


Y 
| 


We combine equations (20) and(21) to form a single bottom boundary 


condition 


—__—_ —_ a SESE 


CZeod2 PLZ, de (22) 


Since we have specified the bottom to be of constant velocity 
and density, a known solution to equation (18) for the bottom region 
is 


Zu Gens: ii coe B® ) (23) 


where Kp, a real constant, is the imaginary bottom wave number, 
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ny ; (24) 


If the amount of energy represented by the bottom solution is to be 
finite, the limit of the integral of Z,, from the bottom boundary to 
infinity must in turn be finite. This condition requires that the expo- 


nential growth term of Zy be suppressed. Thus, (23) becomes 


~Kaz : 
Gage * , (25) 


and (22) becomes 


AaB. & fyy_ 
Zn oz Cb Ve E (26) 


We have now specified two boundary conditions which will 
enable us to find particular solutions to the vertical differential wave 
equation (18). Before proceeding with the technique of finding such 
solutions, let us discuss some of their properties. 

e. ihe Nature of the Solution 

Before we find the solutions, Z(z), to the differential equa- 
tion (18) and boundary conditions (equations 19 and 26), it will be 
beneficial to eunutaes the form and properties of the solutions we 
expect. 

a. In the previous section we saw that the solution in the 
bottom region is of the form of a decaying exponential (equation 25). 


This solution requires that the quantity Kp in equation 25 must be 


real. This requirement in turn implies that 
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E<Vb (27) 


or 


k,>W. (28) 
Cc 
b 


In order to have propagating waves, it must also hold that 


E)0, (29) 
or 


G) 
E 


k,.< , (30) 


min 


Thus, we have an upper and lower bound for E and k, 


V2 £0 (3) 
and 


a Ss eee. 
Cc i Cc 
min b 


(32) 

b. Within these limits the boundary value problem is so con- 
strained that it can‘only be solved in terms of a finite number of 
eigenvalues, Ey and SOS nemckine functions ZA (Z)- This eigen- 
function, properly normalized, forms what is termed a normal mode. 


At sufficient range a vertical sound pressure profile can be approxi- 


mated by a linear combination of such eigenfunctions, 


N 
per = 2, Cm Z mle) . (33) 
m= 


(es) 
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ee wt depths such that E> V(z) the eigenfunction Zn (2) is 


oscillatory in character. We can represent such a function by 


Zz) > Se Sin(Seer) (34) 


where 5 (z) and S(z) are undefined functions of depth. We shall 
refer to such depths as the oscillatory or "real'' region (a reference 
to the fact that the vertical wave number, k is a real number). At 


depths such that Es V(z) the eigenfunction Z,,(z) becomes quasi- 


exponential in character, with a form we can represent by 


Ziay Sx) [Ce2™ « De SJ oe 


where, again, & (z) and S(z) are undefined functions. We shall refer 
to such depths as the exponential or "imaginary" region. 

The depth at which the sign of [E__-V(z)|changes is termed 
a turning point. At the turning point the character of the eigen- 
function changes from oscillatory to quasi-exponential (Fig. 3). 

d. Each mode.corresponds to one or more rays which traces 

paths within the oscillatory region between turning points, surface, 
or bottom boundaries. The angle the ray makes with the horizontal 


is defined by 


Cos hte) = Cte) he (36) 


re) 


From this we see that the mode turning point corresponds to the ray 
vertex depth (that depth at which a ray reaches its maximum depth 


excursion). When the oscillatory region of a mode is bounded by the 


toe 





free surface boundary or the bottom, the corresponding ray is sur- 
face Or bottom reflected, respectively. 

e. The lower boundary placed on k, represented by equation (28) 
has a more familiar and perhaps more satisfying physical meaning. 
Substituting equation (36) into (28) we have 

Cos v) eae 1, 
Cb 
At the bottom this is the expression for the critical angle. Thus, we 
are limiting outselves to consider only those modes whose corres- 
ponding ray strikes the bottom at a grazing angle less than the criti- 
ealangle. Such modes are widely called 'unattenuated modes.'"! 


Also, there exists an infinite set of modes such that 


an (38) 


f 


However, because of bottom reflection loss, those modes tend to 
attenuate rapidly with range. For propagation problems at a con- 
siderable range, ene "attenuated modes'' contribute an insignificant 
amount to the acoustic pressure field, and are ignored. 


f. Given the oscillatory nature of the eigenfunction Z we see 


mM?’ 
that the eigenfunctions for each successive mode must change sign 
one additional time (Figure 4). Each sign change will be referred 


to as a mode crossing. Thus, corresponding to each eigenvalue En? 


there corresponds an eigenfunction 2. with m-1 mode crossings. 
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4, Mode Parameters 
We should now introduce two sets of parameters, derived 
elsewhere, which are used in normal mode calculations. 
a. Phase and Group Speed 
Phase speed describes the horizontal speed of the wave - 
front represented by a mode. As the speed of advance of a wave- 


front of constant phase, the phase speed is given by 


C. eae ( 39) 
Pp <——\* 
Ry 
Group speed is the rate of energy transport in the hori- 
zontal direction. Tolstoy and Clay (1966) give an expression for 
group speed based upon a theorem by Biot (1957). This method cal- 


culates group speed by 


Ca ee (40) 


Oo 


where Vis the normalization integral, 


v= feat de, Sao 


and 


(2) 
o=/e =—— dz, 
-20 C (2>* (42) 


b. Bottom Decay Coefficient 
Kornhauser and Raney (1955) give an expression which 
calculates the effect of a bottom absorption on the mode amplitude. 


Assume the wave number in the bottom is complex and given by 
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iy? Sg (43) 


where B is a bottom attenuation coefficient. Then the horizontal 


wave number must also be complex and given by 
Kr= hy + ia. (44) 


Ibe e is small with respect to & , then the ratio of & to @ for a given 


Cc 
b 
mode is approximated by 
: foe) 
Bee | Bede. as 


me LHE WENTZEL-KRAMERS-BRILLOUIN (WKB) APPROXIMATION 

The WKB approximation is a method borrowed from quantum 
mechanics which will enable us to approximate the eigenvalues, E_, 
of equation (18). In the WKB approximation we assume a form for 
the mode solution which is valid at depths removed from the turning 
points. We then develop a characteristic equation for the eigenvalue, 
based upon the assumed solution. This section follows in many re- 
spects the ae ae given by Tolstoy and Clay (1966) and Lauvstad 
met !). 

1. The Characteristic Equation 

_ Based upon the sign of [E - V(z)] the vertical wave equation 

(18) applies to two domains, which we have termed the "imaginary"! 
region and "'real'' region. Let us rephrase equation (18) into two 


separate equations for the two domains. So doing, we have 
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re he ea, for) Eerie (46) 


and 
22 ee. = 0 for E < Vey. (47) 
Here we have defined 


ee |: ~Vias] when EL? Vw, (48) 


and 
Yo 
K; > [Ves E | when E <¢ \V/«2). (49) 


Now let us assume a form for the solution to equation (46), 
Similar to that which we have previously postulated in equation (34). 


The assumed solution is 


Tewewsieue ~ s. (50) 


By substituting equation (50) into (46) we have as a real part 


Pp p ] 
Sen - Gay) (She wh?) =O, oe 
and an imaginary part 
50a) + 25%2a)S%zy = 0. (52) 


This equation immediately yields an expression for 5 (z), 


Siz) > Af Sz . (53) 


oc 





Now let us simplify equation (51) with an assumption. Assume that 
the first term of equation (51) is negligible with respect to the others, 


or 


Sie 


5 (z) 








With this assurnption we obtain as an expression for S(z) 


Saye wpe ; i) 

or upon integrating 
z 
ena) = t [hyde (56) 
Zo 

We now have as the assumed solution 

Pate: dy,” [a Cos Si (@) + b Sin Sica), (57) 
fr, in a more convenient form, 

Z1(z) = k;* a Sin (Siczy+Qo)] (58) 


Similarly we obtain as a corresponding solution for the 


"imaginary'' region, 


= W-2 S2tz) -S2() 
Zp (2) Ee [Ce De Ih (59) 
where S,(Z) is defined by 


z 
SES -{, dz ; (60) 


Zo 
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Having written the two solutions, let us take a closer look 
at the assumption implied by equation (54). With equations (53) and 


(55) we can rewrite (54) for the "real" region as 


4 | d_ tog hs | «1 (61) 


and for the "imaginary'' region, 





1 d_ tog K 
* SI = iy 


dz (62) 


La 


This assumption requires that the order of magnitude of the vertical 
wave number, and thus the sound speed, vary slowly with respect to 
depth. This is generally the case in the ocean, where the total sound 
speed variation with depth is on the order of five percent. We also 
see that our solutions lose any validity as z approaches a turning 
point. Recalling the analogy of the turning point to the ray vertex 
depth, we see that as a wave approaches the turning point its orienta- 
tion becomes horizontal. Near a turning point the vertical wavenum- 
bers, k, and K,, approach zero and the conditions of equations (61) 
and (62) are no longer satisfied. Thus the solutions represented by 
equations (58) and (59) become singular at the turning point location. 
The WKB method does not, at this point, appear satisfactory 
for the computation of the mode profile or eigenfunction Ze Ae 
ever, if the WKB approximation can be used to obtain an accurate 


estimate of the eigenvalues, E_., a finite difference scheme can then 


Ee used to calculate the eigenfunction, Z_.- 
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Consider a typical V(z) profile, with upper and lower turning 
points, as well as bottom and pressure release boundaries (Figure 5). 
Remember that we have specified a derivative boundary condition at 
the bottom. This boundary condition determines the eigenfunction 
(except for a constant of multiplication) in the ''imaginary''region 
ZT[ Zp. At the lower turning point, zyy,, the "imaginary" region 
solution, Z5(z), corresponds to a particular ''real'!' region solution, 
Z 4 (z). We can designate such a solution by the use of an initial phase 


angle, 97, such that 


(z) ~So(z) 
ks Ce? +De ~* )— A Sin(Sitz)+O@ 63 
where 
ZrL 
So(z) = Soaks : (64) 
and 


S1@) = fohyde P (65) 


Sri 


In a similar manner we can stipulate that in order to satisfy the free 
surface boundary condition, Z(0) = 0, the ''real'' region solution, Zy5 
must correspond to a particular "imaginary'' region solution, Z2. 


Again we can represent this as a phase angle, 0 in the ''real'' region 


S? 
solution, Z,. Thus we have required the argument of Sin (S,(z) +9, ) 


at the upper turning point (or surface) to be our specified value 
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tru 
Shin) elt =f hydz < 0,=O> , (66) 


ZT 
If we let 9,, = - 0, we can rewrite this equation as 
[x dz +0, +0u fe Vee, eee 
Al 
Since this equation describes the arguments for which the sine of the 
integrated vertical wave number is zero (a Dirichlet boundary condi- 
tion), the surface boundary condition will also be satisfied if 
fru 

{ [E- Vel ® de +Q,+ Qu = mar tet2,..,N. (68) 

Zui 
Values of E satisfying this characteristic equation will be the eigen- 


values, E for which we are searching. Note that this equation is 


m’ 
also that given by Tolstoy and Clay (1966) as the characteristic equa- 
tion for stratified acoustic waveguides. 

Now that we have a characteristic equation, some method of 
evaluating 8; and 90,, is required. The WKB approximation provides 
a method for finding the eigenvalues, Ew once 0, and OF are evalu- 


ated. 
@ tlurnine Point Connection 
The most vexing problem in using the WKB approximation is 
connecting the two solutions (58) and (59) across the turning points. A 
number of techniques have been developed, and here we will consider 
two such techniques which we will term the asymptotic and Airy phase 


methods. 
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a. Asymptotic Method 
We have two corresponding solutions represented by 
equations (58) and (59) which correspond to the same function on either 
side ofa turning point. As the turning point is approached let us re- 
quire that the limit of the first derivative divided by the mode value 
be continuous. This is in fact the same as the fluid-fluid boundary 
condition of equation (22). For the purpose of these turning point 


discussions we will designate the turning point depth, z p? equal to 


T 
zero and corresponding to Zo of equations (56) and (60). Further, 


z will be positive towards the ''real'' region and negative towards the 


"imaginary'! region (Figure 6). Thus, from equation (22) we have 


i Za = ej a 
ae at oh ; 7) 


By taking advantage of the approximation involved in equation (61), we 


get by substitution for the above 


lim hy Cot(Sate) +o) = tim Ky Cee? ?-De=*™ 
Z~+Q, z~0- C eee) x De S2(2) (70) 


Taking the limit we have 


b .- Cot (€.) = Bee (71) 
= Oeete 


It should be noted that this approach is simpler than the Airy phase 


technique (described next) and may not be as accurate. It is included, 
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however, since it provides an interesting comparison with the Airy 
phase method in the computer programs. 
b. The Airy Phase Method 

What follows is an outline of a technique developed by 
V. R. Lauvstad (1971) of the SACLANT ASW Research Centre. 

Lauvstad ween asymptotic solutions of the same 
form as equations (57) and (59) in regions remote from the turning 
points. The gradient of the square of the vertical wave number a- 


cross the turning point is assumed to be linear, such that 


[k,] = 02 +... (72) 


This reduces the vertical Helmholtz equation to 


mn + az 7% «QQ, (73) 
One solution to this equation is the Airy phase function 


Z = AA CSe@) + BBi(-Se), ) 


Here 5(z), a depth integration function, is positive in the "real! 
region, and negative in the "imaginary'' region. The asymptotic 
forms of the Airy phase are then compared with our corresponding 
assumed solutions. We can represent the asymptotic form of the 


Airy phase for negative real argument (the "real" region) as 


Zit ee Sizy + (A+B) Cos Scz)| 5 (75) 
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and for positive real argument as 


By matching the above coefficients with those in our solutions we have 
a set of coefficient matching equations representing the required turn- 


ing point connection. 


(77) 
: = AEG 
CL a(2D+C) b As I. 
=1(a-b D=A(a+b) , 
C =A (a-b) aS ) (78) 
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Il. THE 2SROGcekRAM Ma tries 


Here we will discuss the basic methods used by the three pro- 
grams (NORMOI, NORMO4, and NORMOS) presented as a part of 
this thesis. This section is intended as only a brief introduction to 
the programs and their characteristics. A following section will give 


an outline of the program algorithms. 


A. NORMO!I - A FINITE DIFFERENCE TECHNIQUE 

NORMO1 is essentially an adaptation of two shallow water pro- 
grams developed by Kanabis (1972) and Newman and Ingenito (1972). 
For a given value of the horizontal wave number, k;,,or E, the pro- 
gram begins with the bottom boundary condition (equation 22), and 
steps through a finite difference scheme to the surface. The finite 
difference scheme, which is derived in appendix A, is from Kanabis 
(1972). The finite difference scheme is a third-order Newton forward- 
difference scheme which iterates both the eigenfunction and its deriv- 
oie. Using this difference scheme, a search is made for those 
values of k,, or E, which most closely satisfy the surface boundary 
condition (equation 19). Such values are the sought eigenvalues. 

A shortcoming of this technique is that under certain circum- 
stances the eigenfunction at the surface, Z (0): tends to grow expo- 
nentially rather than to decay towards zero. Some mention of this 


phenomenon is made by both Kanabis (1972) and Newman and Igenito 
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(1972). Kanabis recommends decreasing the incremental search 
limits for the eigenvalue Which imposes a consequential increase in 
computation time) to overcome the phenomenon. Newman and Igenito 
set the mode value equal to zero at depths above the last mode cross- 
ing or minimum. It happens that the circumstances for this expo- 
nential growth, or degenerate solution, occur frequently in deep 
ocean acoustic profiles. Thus, a more detailed treatment of the 
phenomenon is required. 

To understand this degenerate solution let us consider a typical 
deep ocean sound profile and our general forms for the eigenfunction 
(equations 34 and 35). Consider a mode whose eigenvalue, ae inter - 
sects the velocity function, V(z), at some depth, Zqy> below the sur- 
face. This is a common (in fact, the usual) occurrence for the lower 
modes inthe deep ocean, where a deep sound channel and thermocline 
usually exist. As noted by Schiff (1955) and Lauvstad (1971), if the 
value of S(z) in equation (34) is not exactly 7 at the turning point, the 

. 4 
coefficients C and D of equation (35) are both non-zero. Thus, the 
exponential growth term has some finite value. Given both sufficient 
depth ie -cun the turning point and the surface, plus sufficiently 
large values of [V(z)-E], the growth exponential will eventually be- 
come dominant. Thus, the solution tends to grow exponentially as the 
surface is approached. 

This degenerate solution can complicate the search for eigen- 


values if it is based solely upon the value of Z(0). With the 
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degenerate solution, the value of Z(0) will oscillate between extreme- 
ly large positive and negative numbers, even with only small changes 
in the horizontal wave number or E. For the lowest modes of a deep 
sound channel profile the values of Z({0) can so oscillate with incre- 


14 Thus, a method such 


mental changes in E of only one part in 10° 
as regula falsi (described later) can become inappropriate since a 
map of Z(0) as a function of E can appear discontinuous (Figure 7). 
In order to circumvent this difficulty, NORMO1 uses a halving 
procedure based upon the number of mode crossings. When the sur- 
face boundary condition is met, the eigenvalue En is the largest 
possible value of E such that the eigenfunction has m-1 mode cross- 
ings. Based on this property, the number of mode crossings is used 
to converge upon the mode eigenvalue. The method of successive 
halving, although not sophisticated, has the versatility to deal suc- 
cesfully with the degenerate solution and any arbitrary profile. Once 
an eigenvalue is found, the function Z, is then checked for a degen- 
erate solution near a surface. If it is present a new profile is 
eee eel in the upper "imaginary'' region and matched to the lower 
profile. This procedure may cause a discontinuity in the derivative 
of Z,, at the turning point; however, its ill effects would usually be 


less serious than those of the degenerate solution. 


B. NORMO4 AND NORMOS5 - THE WKB APPROACH 
In both NORMO4 and NORMOS the vertical wave number is inte- 


grated between the upper turning point (or surface) and lower turning 
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point (or bottom). The programs differ only in the turning point 
connection formulae. Both programs integrate the vertical wave 
number using the trapezoidal rule over the depth grid points. The 
eigenvalue E_ is converged upon using the regula falsi method. 
This method makes successive estimates of the zero value of a func- 
tion by the linear interpolation of a negative and positive value. The 
method is represented by 

Sy) Ee eee AS}... (79) 

(AS. - AS-) 

where the subscripts + and - correspond to the values associated 
with the last positive and negative values of MS; and AS, the differ- 


ence function, is defined by 
Zev 
AS = {hy dz +8, °*0. -ma - (80) 
zn: 


If the mode is not found after a set number of iterations, the regula 
falsi method is dropped for a halving process. 

The phase effect of the upper and lower turning points, 9,, and 
Or» D6 evaluated using the asymptotic method for NORMO4 and the 
Airy phase method for NORMO5. The vertical wave number is integ- 
rated from the uppermost turning point (or surface, if surface re- 
flected) to the lowest turning point (or bottom, if bottom reflected). 
Both programs must then provide, in addition to the values of @,, and 
Or» the effect of any intermediate "imaginary" region or barrier. 


This occurs when multiple sound channels are encountered (Figure 8). 
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Details of the derivation of the specific formulae for oF 6,, and the 
phase effect of an intermediate barrier are given in appendices B and 
C for NORMO4 and NORMO5 respectively. 

If the mode is surface reflected or bottom reflected, both pro- 
grams use the same formulae for the phase effect of the reflection, 
0, and @;. Since the free surface requires Z(0) = 0, for the surface 


reflected wave, we have from equations (63) and (66) 
fru 
J bade +0, =8.= mr (8.1) 


and thus 


Deo. (82) 


For the bottom reflected mode, we require the conditions ofa fluid- 
fluid boundary to be met (equation 22). Substituting equations (63) 


and (25) into (22) we have 


hz Cot @.)= - & Ke , (83) 
or 


(84) 


=Arcla Cb de 
6. A ct n (Se ), 


where k, is evaluated at the water side of the bottom boundary. 
The occurrence of multiple sound channels or ducts makes the 
WKB method more complicated. If the intermediate "imaginary" 


region or barrier is of sufficient extent and magnitude, sound 
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propagation in each channel becomes independent of the other. When 
this is the case, the waves in one channel undergo an effectively total 
reflection by the barrier. Thus, complete sound ducting occurs 
entirely within that channel. Lauvstad (1971) gives the ratio of the 


amplitudes of the reflected and incident waves ina duct as 


R,= e -et . (85) 


eo ae oe 
where L is the integral of the imaginary vertical wave number across 


the barrier, 
- (86) 
| = [Bade ° 
ot 


This result can be derived using either the asymptotic or Airy phase 
connection formulae. As L increases, R, rapidly approaches unity, . 
and near total reflection is achieved by the barrier. The WKB pro- 
grams NORM0O4 and NORMO5 evaluate L, and if it is of sufficient 
magnitude, the programs treat the mode propagation in each channel 
independently. The programs accommodate a maximum of one main 
sound channel (which contains the velocity minimum), and one upper 


and one lower sound channel. 
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IV. THE BASIC ALGORITHMS 


A. NORMO1 
Po Geteral Description 

NORMO1 iterates the wave equation vertically from the 
bottom to the surface using a finite difference scheme over a depth 
grid of up to 501 points. The program is written in IBM FORTRAN 
IV with subroutines used as major functional blocks. This modular 
programming is designed for both ease of modification and debugging. 
A description of the sequence of events follows. 

NORMO!1 is basically composed of three nested loops: a 
frequency loop, a mode loop, and an iteration loop. The input data 
are read using an input ene. The required input data are 
listed in appendix D. Once the data are read and appropriate con- 
stants computed, the sound velocity profile is linearly interpolated 
for all grid points. Here the frequency loop begins, the appropriate 
frequency is selected, and the velocity function V(z) is computed for 
the grid points. Next the maximum and minimum horizontal wave 
numbers and maximum number of modes are computed. If the fre- 
quency is below cut-off (the maximum number of modes is zero), a 
new frequency is selected. At this point the mode loop begins, and 
the following steps are repeated for the first to the highest requested 


or possible mode. 
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First, a method of successive halving is used to find two 
values of the horizontal wave number k, which have m-l and m mode 
crossings. The desired eigenvalue must lie between these two values. 
Then, again using a halving process, the program converges upon the 
eigenvalue k.. The mode is considered found when one of two con- 


ditions (regulated by an input parameter, IEX) is satisfied; 


=IEX 
AC Ate (0, (87) 


or the horizontal wave number changes between iterations by an 


amount such that 


Ree, a he 
Re, 


4 


4 








Once the mode is found, the eigenfunction is checked for a 
degenerate solution (exponential growth near the surface). Ifa de- 
generate solution is present, a new function is iterated from the sur- 
fae to the upper turning point, then joined with the lower function. 
The mode is then normalized to its absolute maximum value, and the 
profile within the homogenous bottom is computed. The computed 
mode is plotted using a printer plot routine. A subroutine then com- 
putes the normalization integral represented by equation (41), and 
the group speed, phase speed and mode attenuation coefficient. At 
this point the mode and frequency loops end. 

Included in NORMO1 is an output subroutine which provides 


a punched card format listing of the eigenvalues, mode parameters, 


=| 








sound velocity profile and mode profile. The output subroutine was 
written to provide local input to CALCOMP plots and propagation 
loss programs. 
2. The Subroutines 
The following is a brief description of the major subroutines 


which comprise the NORMO1 program. 


INPUT 1 provides the input data for the program. A user 
supplied page heading is read on the first card, followed by the run 
parameters, list of frequencies desired, and a maximum of 50 depth 
and sound velocity pairs. If the parameter representing the number 
of depth and velocity pairs, NUMV, is negative, all length measure- 
feents read are converted from feet to meters. This subroutine also 
prints the input data for reference. 

INTRPL linearly interpolates the input depth and sound 
velocity pairs to compute the grid point values of sound velocity. If 
required, this subroutine also extrapolates the last sound velocity 
to the bottom assuming an isothermal, isohaline gradient (0. 017sec~. 

MAXMIN computes and checks the maximum and minimum 
values of the horizontal wave number, a finds the number of modes 
represented by the minimum horizontal wave number (the maximum 
number of modes). If the maximum number of modes is zero, the 
frequency is below cut-off and the program selects a new frequency 
or terminates the calculations. If the maximum possible number of 


modes is below that requested, this parameteris changed and noted. 
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HALF searches for two values of the horizontal wave number 
with m (the mode being searched) and m-1 mode crossings. These 
values represent lower and upper bounds on k, which are then suc- 
cessively narrowed. At the same time this subroutine 'maps'" the 
values of k, versus mode number in order to facilitate the search 
for later modes. 

ITERAT, the heart of the program, iterates the eigenfunc- 
tion from the bottom to the surface using equations (97) and (101). 

In addition, this subroutine counts the number of mode crossings and 
finds the maximum absolute value of the eigenfunction. 

DNORM1 normalizes the eigenfunction with respect to its 
maximum absolute value. 

BOTTOM creates an exponentially decaying bottom mode 
profile from the bottom to an input depth, PLTMAX. 

MODPLT and CZPLOT are printer plot subroutines which 
plot the mode shape and sound velocity profile. Both these subrou- 
tines use a Naval Postgraduate School utility printer plot subroutine, 
UTPLOT. 

FIXUP tests for a degenerate solution near the surface. If 
the degenerate solution exists, this subroutine applies the finite dif- 
ference iteration scheme starting from the surface downward to the 
upper turning point. This upper mode shape is then scaled to join 
the original lower function at the turning point. If required, a new 


value for the absolute maximum of the eigenfunction is then found. 
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INTEGR computes the mode normalization integral (equa- 
tion 41) and the integral represented by (equation 42). The phase 
velocity is computed; then the two integrals are used to compute the 
group velocity. Finally, the mode attenuation coefficient (equation 


45) is calculated. 


B. NORMO4 AND NORMO5 

1. General Description 

NORMO4 and NORMO5 are identical with the exception of the 

turning point connection formulae (appendices B and C explain these 
differences). Similar to NORMOI, they consist of three nested loops; 
a frequency loop, a mode loop, and an iterative loop represented by 
the subroutine SEARCH. In addition, the input, intoeeotetion, plot - 
ting, maximum/minimum and mode integration functions are per- 
formed by subroutines similar to those used in NORMOL. 

The sequence of events in the WKB method programs is 
similar Z NORMO1. NORMO4 and NORMO5 read the input data and 
then interpolate the sound velocity profile over a maximum of 1001 
ities, grid points. If requested, a sound velocity printer plot or 
punched card format sound profile is provided. The frequency loop 
then computes the velocity function V(z) and checks for maximum 
horizontal wave number, minimum Hone ottesl wave number, fre- 
quency cut-off, and the maximum number of modes available. At 
this point the mode loop begins and the search for the eigenvalue is 


controlled by the SEARCH subroutine. 
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After the mode is found, the mode shape is computed using 
a finite difference scheme. The finite difference scheme is fitted be- 
tween the upper and lower turning points, then a solution which decays 
away from the turning points is calculated. After the mode shape has 
been calculated and normalized, the mode is integrated to provide the 
normalization integral, phase and group speed, and the mode attenu- 
ation coefficient. Finally, if requested, a printer plot of the mode 
shape is provided. At this point the mode and frequency loops end. 

Because of their importance to these programs, the SEARCH 
and WKB subroutines deserve more extensive description. 

2. SEARCH Subroutine 
The SEARCH subroutine controls the iterative search for the 
elgenvalues kK The subroutine first makes an initial guess of an 
upper and lower bound for the eigenvalue. The WKB subroutine is 
then called to compute the difference function given by equation (80). 
If the difference function for the lower bound is not positive, the 
bound is successively relaxed until a positive value is found. 

The next guess for the eigenvalue is calculated using a 
regula falsi Se technique ras equation 79). The re- 
sulting incremental change in eigenvalue and the absolute value of the 
difference function is checked to test whether the desired accuracy 
has been achieved. The WKB subroutine is then called to find the 
difference function for the new trial eigenvalue. If the resulting dif- 


ference function is negative the upper bound is replaced; if it is 


oe) 





positive the lower bound is replaced. A new guess for k, is then 
estimated, and this iterative process continues until one of two con- 
ditions is met. The mode is considered found when the difference 


function falls below an absolute limit, 


ZS leatioe ** © (89) 


or the value of the horizontal wave number approaches machine 


accuracy, 


Atl a 


he 
a 


Ve jee (90) 








If after a set number of iterations, the regula falsi method has failed 
to converge upon a solution, it is replaced by a cruder and slower 
(but more versatile) halving process. 

If more than one sound channel exists, the program must 
determine if the two channels are independent of one another. If the 
L of equation (86) is of sufficient magnitude, the subroutine SEARCH 
computes the eigenvalues for each channel separately. Thus, after 
an eigenvalue has been found, two tests are made. First, it is deter- 
mined whether an upper or lower (separated) sound channel exists. 
Then, if an upper or lower channel exists, it is tested to determine 
whether it will propagate an additional mode. If such a mode can be 
propagated, the SEARCH subroutine computes its next mode within 


that channel. 
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3. WKB Subroutine 

The purpose of the WKB subroutine is to integrate the verti- 
cal wave number between the upper and lower turning points. The 
subroutine consists of a loop which steps through the depth grid from 
the surface to the bottom. Inany "imaginary" region the imaginary 
vertical wave number K, is integrated for use with the turning point 
connection formulae. Each time the depth loop exits an "imaginary"! 
mecven for a 'real'' region, either the upper turning point phase 8. 
or the phase effect of a barrier must be calculated. If the barrier is 
considered to be separated (if L in equation 86 is larger than a set 
value) the integration is stopped with the upper channel. 

When the depth loop reaches the bottom, the bottom phase, 
Or: is computed for either a bottom reflected or a decaying solution 
as appropriate. Finally, the difference function is calculated and 


program control returns to the calling program. 
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Ve RESU LIS 


All three programs were run with a variety of sound velocity 
profiles and the calculated eigenvalues examined. The profiles select- 
ed include the published results of other models and profiles with 
known analytic solutions. In all the test runs the input variable IEX 
was setat 10. Thus, the iterations were halted when the absolute 
mewative value of the eigenfunction at the surface, Z(g9) , or the dif- 
ference function, A S , (equation 80) achieved a value less than tO meee 
Additionally, computation was halted if the eigenvalue was incremented 


14 


less than 10 ~* during a halving process (in other words, when the 


solution approached machine accuracy). 


Pee lrkol CASE ONE - NEGATIVE GRADIENT 

This test case is one given by Kanabis (1972). The profile is 
characterized by a simple negative velocity gradient ina shallow 
water channel (Figure 9). The profile has a fast fluid bottom of the 
same density as the water. The test results are computed for a fre- 
quency of 1000 Hz. Table 1 lists the eigenvalues given by Kanabis 
and computed by NORMO1, NORMO4 and NORMO5 in ae The list 
includes a set computed by the Kanabis program and a set computed 
by an unpublished program by C. L. Bartberger and L. L. Ackler. 
Except for the higher modes of NORMO5, the results from all three 


programs agree to the fourth decimal place or one part in 10°. 
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Table 2 gives the NORMOI1, NORMO4 and NORMOS5 results to 
the accuracy computed and relative differences between NORMOI1 and 
the other two programs. NORMOl1 required 56.23 sec of central 
processing unit (CPU) time and 412 iterations to arrive at the solu- 
tions. NORMO4 and NORMOSd required 23.50 seconds with 83 itera- 
tions and 22.74 seconds with 79 iterations respectively. The CPU 
times cited include the time required for all operations including in- 


put, sound speed profile manipulation and graph printing. 


me Laks T CASE TWO - POSITIVE GRADIENT 

The second test case is a shallow water example cited by Newman 
and Ingenito (1972). The profile is characterized by a slow, variable 
positive sound speed gradient over a fast, dense bottom (Figure 11). 
The computation frequency was 700 Hz. It should be noted that the 
accuracy parameter used by Newman and Ingenito was an absolute 
value of less than 0.001 for the surface eigenfunction Z(0). This 
roughly corresponds to a value of IEX between 3 and 4 in the NORMO 
programs. With the exception of NORMO5, all the programs agreed 
to at least the third decimal place (Table 3). Considering the accu- 
racy limit set for the Newman and Ingenito run, this agreement is 
about as good as may be expected. NORMO1 required 27.43 seconds 
and 196 iterations to complete computations, while NORMO4 and 
NORMOS required 14.73 seconds for 77 iterations and 14.92 seconds 


for 79 iterations, respectively. 
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fe LEST CASE THREE SoeMMERIC ee buUCr 
A profile was constructed corresponding to a symmetric wave 


duct where the sound speed is given by 


4 2k _ 2 


between two depths. The solution to this profile is given by Tolstoy 
and Clay (1966) as 

,,, = [2 - ora. (@m-L}] (116) 

m Ce 

The actual profile dimensions, shown in figure 12, are taken from 
Tolstoy and Clay. The profile was tested at 30.0 and 60.0 hertz. 
The computed eigenvalues agree with the anlytic solutions to one 
part in 10°, with the exception of the higher modes of NORMOI1. 
Table 4 shows comparisons between the analytic solution and 
NORMOI!, NORMO4 or NORMO5 computed values. 

Of significant interest is the fact that the errors appeared rela- 
tively stable for NORMO4 and NORMOS5. The errors for NORMO4 
are consistently 0.45x107° to: 0. Are ie meters -l lowfor 30 Hz, and 
0.10x107° “ts 0. 12x107> meters — low for 60 Hz. In contrast, the 
errors for NORMO1 increased from 0. 39x107° meters ~! low to 
6 


12.5x10~° meters 7! high for 30 Hz, and from 0. g0x10~° to Zao On 


meters ~! low for 60 Hz. Figure 14 displays a graph of the program 


errors. Figure 13 shows the first three modes. 
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m. TEST CASE FOUR - DEEP SOUND CHANNEL 

Test case four is a typical deep ocean sound channel. The pro- 
file is characterized by a deep and sharp sound channel and by a 
bottom of the same density as the water (Figure 15). With the excep- 
tion of the first mode, the eigenvalues agree within one part in i 
(Table 5). Note that for the higher frequency, 60 Hz, the agreement 
is somewhat improved. It is interesting that the first mode should 
show a larger disagreement among the programs. This effect may 
be a result of the combination of the sharp sound channel axis anda 


grid spacing of 8 meters. 
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Frequency: 30 hertz 
MODE NORMO1 
] 0. 12614128 
(i S18) if Hols 
3 53 (a6 6 
4 514437 
5 492671 

6 472443 
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9 419601 
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Frequency: 60 hertz 
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-32.84 0.12610844 


74 


DIFF 
NORMO4~ x10°® 

0. 12610844 
567933 - 1.97 
538187 8.00 
513923 -5.14 
Ae 7 4 AAT 
ZL (7D Cn || Ze 
453659 + - 0.86 
436012 - 1.50 
419117 + - 4.84 
400638 -34.68 

~45, 35 

196827 1. 68 
154783 1.25 
(118986 7p yl 
089755 5.78: 
063571 - 0. 83 
039343 -2.18 
016618 0.24 

0.24994890 iG aS 
974255 Seal 

TABLE V 
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mee LEST CASE FIVE - DOUBLE CHANNEL 

Test case five is a double sound channel of nearly symmetric 
dimensions (Figure 16). This profile was used to test the multi-duct 
properties of the NORMO4 and NORMOS5 programs. Note that for 
30 Hz thedifferences between NORMO1]1 and NORMO4 (or NORMO5) 
increase by almost an order of magnitude for modes 5 and 7 (Table 6). 
For these modes the barrier connection formula (equation 114) is 
employed and the ducts are considered "connected.'' The mode pro- 
files (Figure 17) for these modes also do not completely correspond. 
Modes 5 and 7 appear to be "upside-down" in NORMO4, with the 
dominant amplitude waves in the lower channel. This is due to the 
fact that NORMO4 and NORMOS5 iterate the vertical wave number 
from the surface downward. Because of this, the phase at the upper 
turning point of the lower channel represents the phase angle for the 
third and fourth mode in the lower channel. Equation (114) remains 


correct, although here it has been incorrectly applied. 
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RESULTSJOF TES) CAGE stick 
DIFF 
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MODE NORMOI! NORMO4 
frequency: 30 hertz 
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Vie CONCLUSIONS 


A. GENERAL 

As seen by table 7, the WKB method has performed consistently 
faster than the iterative technique. The WKB programs (NORM04 
and NORMOS5) found solutions with about one-half the Central Pro- 
cessing Unit (CPU) time required by the finite difference program 
(NORMO1). In addition, one should consider the method used for in- 
tegration of the vertical wave number in the WKB programs. The 
trapezoid rule was used which required the evaluation of a square 
root at each grid point. Although simple, this is a time consuming 
task, while it is not necessarily very accurate. Ifa faster integra- 
ting procedure were substituted for the trapezoid rule, the time re- 
quired for each iteration could be reduced considerably. Thus, in 
addition to showing an improvement in computation time, the WKB 
method also offers the potential of greater computational speed 
through the use of sophisticated integrators. 

The WKB and finite difference methods appear to have similar 
accuracy in cOmparison with other programs. The accuracy attain- 
able is adequate to calculate useful modes for a propagation loss pro- 
gram. However, the errors involved in the differences between 
successive eigenvalues will not permit accurate coherent phasing of 
the modes beyond the first or second convergence zone . Coherent 


phasing at extremely long ranges is, in fact, attainable only with an 


“y, 





CPU TIMES FOR VARIOUS RUNS 


NORMO! NORMO4 NORMO5 
Test Case One Soe 2o sce 23.50 sec 22. (4ecee 
Test Case Two 27.43 12 as 14. 92 
Test Case Three 106. 35 38.42 Soe 
Test Case Four 105.67 40.23 40.43 
Test Case Five Mma, O7 44,53 47.09 
ZOLATL 409.75 sec lol. 4) sec 164.65 sec 
TABLE VII 
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enalytic profile and solution — an ideal mathematical model. (in 
fact, when one considers the accuracy of the available sound velocity 
data, it is difficult to imagine any model attempting to yield accurate 
phasing effects at any large distance. Since the position of a target 
is unknown, in practice it is more important that the amplitude and 
nature of phasing effects be know, than for the sharp interference 
peaks to be accurately predicted. Thus, accuracies of one part in 


10° 


should enable a model both to take account of the phasing effects 
in the first convergence zone, and to give a general description of the 
effects of phase interference at longer ranges. 

Although the WKB programs do not correctly interpret the barrier 
effects, the problem appears to lie in the approach taken in the pro- 
gramming, rather than in the mathematical and physical relationships. 
Within any duct, the effects of any upper and lower ducts must be 
treated as a turning point phase effect similar to the top and bottom 
boundary conditions. The modes 5 and 7 of the NORMO4 30 Hz double 
channel case are unwitting examples of this. In NORMO4 and 
NORMOS the "connected" ducts are considered as a single wave sys- 
tem. In fact, with any sized barrier, no matter how small, one must 
consider the wave system within each duct separately. In this case 
the turning point phase angle between the duct and barrier is analo- 
gous to the impedence due to the other duct felt by the wave system in 


its duct. 
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Boe FUTURE PROGRAM 

It appears that a faster WKB method normal mode program can 
be written which will properly interpret a barrier. In order to pro- 
perly consider the multiple channel cases, a mapping procedure 
should be used before actual mode searches commence. The sound 
profile would be divided into "levels" (Figure 18), found by selecting 
all the velocity maxima and minima. Each "level'' represents a 
range of the possible eigenvalues, En’ Within each "'level'' there is 
a single combination of ducts, barriers and boundaries. At the start 
of computations the program would determine the limits of each level. 
Then, at the beginning of the frequency loop, the program would de- 
termine the maximum and minimum number of modes for each duct 
within each level. Thus, equipped with a map of the profile charac- 
teristics, the program would compute eigenvalues for each duct and 
level in sequence (figure 19). After all the eigenvalues are computed 
and sorted (some may be out of numerical order), the eigenfunctions 
can be computed quickly using the fimmte difference procedure of 
NORMO1. Geebined with a faster integrating method, sucha pro- 
gram should yield fast, accurate normal modes for an arbitrary 


sound speed profile. 
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ped eal NIB PCy" 


NORMO!] FINITE DIFFERENCE SCHEME 


In order to iterate equation (18) through the depth grid, a finite 
difference scheme is used which computes values of the eigenfunction, 
Z(z) and its derivative, Z'(z) at each grid point. This scheme was 
originally developed by the Arthur D. Little Company [Report No. C- 
70673, 31 January 1969]. 


Z(z) is expanded into a Taylor series, 


Bierh) = Zz) + hAta + i “z) oa > (91) 
where h is the depth grid spacing. From equation (18) we have 


aX ays - [E- Vea)] Biz) = — Fezy4cz) (22) 


where F(z) represents the square of the vertical wave number. Using 


a forward difference for Z'''(z), 


ay) = Aah) (93) 
: h 


Substituting equation (92), we have 


ees) > 4 [FnZee = Fiz+h)Bczeh)| F ve 


Combining equations (91), (92) and (94), 


Zieehy® Berrh E%a)- sf Fens + We lFin Zen Fee Bees] , (95) 


28 





we have on clearing, 
Z zsh) E : IF Fash) : Zce)|t- he Fea) Ba ual (96) 
3 


Placing the above result in a vector subscript notation we have the 


formula used to compute the succeeding eigenfunction, 


Bien * Bilt 9 hu) 2% a 


(4 PFs) 
We must now compute the value of the derivative of the eigenfunction 
Z'(z) at the next grid point. In order to do this let us express the 
eigenfunction at our original grid point as a backwards Taylor series 


of the new grid point, zth, 
Za) = Zlerh) -hB (e+h) ae Z. (eh) - ~ Dente. (98) 
Substituting equations (92) and (93), we have 
Ate) > Zizrh)-h2’a+h) “ Fearhy Zteeh) ~ LiFe 2 Fea Zeerh) » (99) 


Rearranging, we have, 


2 ; > 
Mech) = 4| Zeaem ilole Fcx-h)} -Zas(1+% Fe) ’ (100) 


or in vector notation, 


fi = ( = et) = Zoe + Ht re 
Bir = 4 [Biss Gi 3 Fi ) Zz (1 YF) | ; (101) 
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APPENDIX B 


NORMO4 CONNECTION FORMULAE 


NORMO4 requires the evaluation of the phase at the upper and 
lower turning points in addition to the phase effect of a barrier. For 
the purposes of this discussion we will align the origin of the depth 
axis with the turning point, designating the positive direction towards 
the ''real'' region. 

Sernsider the upper turning point case (Figure >). The free sur— 
face boundary condition requires that the "imaginary" region solution 


disappear at the surface. From equation (59) this requires that 


Z(-H) =K,CH)| Ces + ee eo 510). (102) 


Thus, we can define C and D (except for a constant of multiplication); 


Cine e ~ S2(-H) ea Sa (-H) (103) 
Applying equation (71), we have 
-~S2(H) _ e S2CH) 
bo uCH , eGR) ” p NES, 
or 
0 
Bu = Arctan [Tanh [Kz de) | . (105) 
-A 








Now, consider the bottom boundary condition and lower turning 
point (Figure 5). The bottom boundary condition (equation 26) re - 


Gu@ires 


fe, Xe 2 Ke Ce S264) _ De s204) (106) 
Cb Ce SeC4) 4 De -seCH) ? 


where Eee is defined as the imaginary vertical wave number evaluated 


at the water side of the bottom interface, 


2 _ a 4 ae 7 
Ke -[Veu) E| 5 ae . (107) 





Thus, the imaginary solution can be written by defining C and D such 


that 


Cathie kee 


108 
D= (eo Ks - @K,)e %o™ a 


Now, by applying the turning point connection formulae (equation 71) 


we have 
1 _(CoKe+ CoN geo + (CbKse- Kaen) 
ey (Chie eee (rns tenner (109) 
or 
8,=Arctan {D2 +4] , (110) 
Ls i 1 
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where D, is defined as 


is = Cora ns Soles 2S2(-H) 
ee a eli 
Coins, + Coie (111) 


Now, consider the situation with an upper and lower channel 
separated by a barrier (Figure 8). Let the z origin be at the lower 


connection point, z,, with the positive direction downward. From 


2 
equations (59), (69) and (71) we have for a connection at the upper 


point, Zy> 


ay > De * + Cet : 
“LL wu 


where L is as defined in equation (86). Letting A, equal the numer- 


ator and By the denominator, we have, except for a constant, 


Gee (a, -b,)e> 
(113) 
Di=5 (a, +b, Je" . 
Applying these expressions to equation (71), we have 
G2 (ay+bijeh + G@i-bje™ (114) 
Ss - (a, +b,eé - (ar-bi)e* 
or 
Ome =eferawern | isn Q:+1)et+(tan@i-t)e*] - (115) 
(tan Ox +1) e> - (tan@:1-i)e 
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APPENDIX C 


NORMOS5 CONNECTION FORMULAE 


For NORMO5 the development of the upper and lower "imaginary'' 


region particular solutions by the application of the free surface and 


bottom boundary conditions is identical to that for NORMO4. The 


differences in the two programs lie in the application of the Airy 


phase connection formulae vice the asymptotic connection formulae. 


To obtain the upper turning point phase, we apply equation (103) to 


equation (77) 


a seals @s2CH) es eae 


ce + @ 82 CH)) ; 


or 
enainrcien oy 
where D, is defined as 


D, ay e 2S2CH) 


Note that as H approaches zero (the turning point 


face) the value of 0, in equation (119) approaches 


(116) 


(7) 


(118) 


(119) 


approaches the sur- 


arctan (1/3) vice the 


phase we required for the free surface boundary condition, (@,, = 0 ). 
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To derive the lower turning point phase, we apply the particular 
lower "imaginary'' region solution (equation 108) to the Airy phase 


connection formulae (equation 77), 


ee | - Yo S CH) | + Xo oheGle 
a - A [2(eKs Go Kg) e7*' “+ (Po Kg + &Kg)e | (120) 
: b= & [2 CoKs- Co Ka) eo -(e, Kes CoH g)e | (121) 
Dividing, we have 
2D2+ 1 
: Eber 1 
6, Ae [e2e 7" (122) 


where D, is defined as in equation (111). Note that this equation also 
does not approach the bottom reflected boundary condition (equation 
84) as H approaches zero (the lower turning point approaches the 
bottom). 

Lauvstad (1971) derived an expression for the phase coefficients 
of a wave entering a barrier region in terms of the phase of the wave 


leaving the opposite side of the barrier: 
O.=2 c +b)e "+ (az - b:)e"| 


(123) 
b, = A [(orsbe)e™ = (on ~be)e*| 


By rearranging, we obtain 


el 





ie: -E[ra- bi) 6 "eee Caa bien] 


(124) 
which agrees with equation (114), the NORMO4 result, 
0. = Arctan (CE es + bi) es (ar- bi)e™ (114) 
(aq ba) eo (a,~- ber" : 


Re 





APPENDIX D 


INEUL BATA 


The following format is used for the input data to all three programs, 


NORMO1, NORMO4 and NORMOS. 


VARIABLE 
CARD ONE 
HED (20) 

CARD TWO 


NUMV 


NMOD 

IE Xx 

N 

IOUT 

ete V * 

CARD THREE 
lol 

CB 

DRO 

DRB 


ee OT 


*NORMO4 and NORMO5 only. 


FORMAT 


20A4 


110 


110 
110 
110 
110 


I2 


P10. 3 
FLO..3 
On 
FIO;3 


Flo..3 


MEANING 


HEADING. TITLE CARD 


RANGE 


Number of Depth/Sound Speed 2 - 50 


Pairs. Negative if conver- 
sion from feet to meters de- ° 


sired. 


Number of Modes Desired 
Iteration limit. 

Number of grid points desired 
Type of output desired 


Output device for card file 


Bottom Depth 


Bottom Sound Speed 


Water density 


Bottom density 


Maximum depth to be plotted 


ve 


1 - 100 
1 - 14 
100 -500 


100-1000* 


see below 


HEC 





CARD FOUR 


NFREQ Ie Number of frequencies 
desired 1 - 20 


@eRD FIVE (SIX) 

meeO (NFREQ) 8Flo.3 Frequencies desired 

CARDS SIX TO 5tNUMV (SEVEN to 6+NUMV) 

DEP (1) eNO Depth DEE i> 
Baier ab 1) 

CC (I) FO. 3 Sound Speed 

For NORMOI, IOUT=!I gives a punched deck output of the sound 


profile and mode parameters and profile. For NORMO4 and 


NORMO5, the IOUT parameter gives the following types of output: 


IOUT PRINTER GRAPE VP UNCHED CARDS ~™ bisk CULrEeUL 
0 YES no no 
] YES YES no 
Z YES no YES 
3 ne | no YES 
4 no YES no 
2) no no no 
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NORMO1! FLOWCHART 


a PRs 
READ AND 
CHECK INPU 


DATA 





V/ 
LINEARLY INTERPOLATE 
SOUND PRO FILE 


SOUND VELOC 
LILY GaeaPH 













\/ 
PUNCH 

\ SOUND VELOC 

ITY PROFILE 
2K 
NS 

— FREQUENCY LOOP BEGINS IFREQ=1, NFREQ 
3 &, 
COMPUTE ViIboOcrl Y 
FUNCTION 

\/ 







CHECK FOR 
MAXIMUM AND 
MINIMUM NUM 
BER OF MODES 








below cut-off freq 


a5 





y 


\/ 


P< SMODE LOOP BEGINS M= 1, NMOD 


USE METHOD OF HALVES 
TO GET WEEE R AND BOW 


ER BOUND OF EIGENVALUE 
WITH M-1 AND M MODE 
CROSSINGS 





\/ 
ITERATION LOOP BEGINS 


4 7 
— COMBUTE NEXT GUESS 
BY HALVING 


MODE LOOP 


\/ 


FREQUENCY LOOP 







INCREMENTAL 
CHANGE < DESIRED 
ACCURACY ? 


yes 


no 


=) Wy, 


CALL ITERAT 


\4 
ITERAT IO 


4 COUNTER 
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SURFACE MODE 
VALUE << ACCURACY 
DESIRED? 


Yes 


yes no 


CROSSINGS < M2? 


REPLACE UPPER 
BOUNDS 





REPLACE LOWER 
BOUNDS 






\4 


END OF ITERATION LOOP 


MODE LOOP 


CHECK FOR DEGENERATE 


SOLUTION AND FIXUP 
IF REQUIRED 


FREQUENCY LOOP 





MODE HAS BEEN FOUND 


j 


NORMALIZE MODE TO ONE 





a7 





FREQUECNY LOOP 


MODE LOOP 


\Z 
COMPUTE A BOTTOM SOLUTION f 


\/ 
GRAPH THE MODE 


t 


COMPUTE NORMALIZATION 
INTEGRAL, GROUP AND 


PHASE SPEEDS 








\/ 
~— PUNCH MODE 
OUTPUT 7 bs OUTPUT DECK 
12 Uy, 
END OF MODE LOOP 
15 “Wy, 
END OF FREQUENCY LOOP 
\/ 





NORMO4 AND NORMOS5 FLOWCHART 


READ AND CHECK 
INPUT DATA 


LINEARLY INTERPOLATE 
THE SOUND VELOCITY 


GRAPH SOUND 
PROFILE | 





PUNCH SOUND 


PROFILE DECK 





no 
—_| 


\) 


— FREQUENCY LOOP IFREQ = 1, NFREQ 
dj Y’ 


a] 


a 





——— 


FREQUENCY LOOP 
- MODE LOOP —— 


—— 


COMPUTE VELOCITY 
FUNCTION V(Z) 


CHECK FOR MAXIMUM 

AND MINIMUM MODES 

below cut-off 
frequency 





\ 


> MODE LOOP MODE = 1, NMODF 


\ 


CALL SEARCH 


(SEARCH FOR 
EIGENVALUE ) 











CALL SHAPE4 


(COMPUTE THE 
MODE SHAPE) 








INTEGRATE MODE, COMPUTE 
PHASE AND GROUP SPEED, 
AND DECAY COEF. 


100 





FREQUENCY LOOP 
MODE LOOP 


NO ves 


y 


PRINT MODE GRAPH MODE 
PARAMETERS 


PUNCH MODE 


SHAPE DECK 





4 \) 
MODE LOOP ENDS 
5 : \) 
| FREQUENCY LOOP ENDS 


J 
Co swe 


101 





SEARCH SUBROUTINE FLOWCHART 


no 


SET SOME 






oe > CONSTANTS _ 


1 
no 
3 J 
— MAKE 
FIRST GUESS 
ICALL WKB fle" 
DUCT 


CUPPER DUCT? 
no 


UPPER DUCT 
HAVE ANOTHER 


Ay v 


set FOR 
SECOND 













<a 


102 


112 \ 


SET UPEER 
LIMIG 














LIMIT 
IN LOWER 
DUCT? 


DIFFERENCE 
»0? 


no 






yes 





y yf YY 


Sal LIMITS 
FOR 
UPPER WELL 


32 









LOWER 
DIFFERENCE 
202 






yes 











no 






STILL 
IN LOWER 
DUET 7 


RELAX 
LOWER 






yes 









no 
5s VJ 
UPPER 
[Eps T ES 6 GO TO 
DUCT? TOP DUCT 







y, 


yes 
10 





UPPER : 
DIFFERENCE 
C50? 






16 \ 










COMPUTE NEXT GUESS 
BY REGULA FALST 








CALL WKB 


7 \ 


ITERATION 
COUNTER 
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en y 










RE PLACE 
UPPER 
LIMITS 












LOWER 
LIMITS 








SIFFERENCE 
< ACCURACY 
DESIRED? 


24 









FIND NEXT 
GUESS BY THE EIGENVALUE 
HAS 


BEEN FOUND 


no 


HALVING 





ITERATIONS 
< ITMAX? 


ho 


“is \ 









UPPER OR 


LOWER DUCT? 
lower 


26 \) 27 \/ 


NEXT DUCT 
WILL BE CENTER 


: 
RECORD LIMITS 









NEXT DUCT 
WILL BE TOP 


RECORD 


LIMITS OF 
CENTER 


SECU eS DUCT 
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| 









1S THERE 
A LOWER DUCT? 






yes 






no 


wt 


CALL WKB V 


> 


NEXT DUCT 
WILL BE TOP 





\ 


WILL LOWER sss \/ 


DUCT TAKE ANOTHER : 
29 y 
NEXT MODE WILL 

BE TOP DUCT 


30S 


SET LOWER 
LIMITS FOR 
LOWER DUCT 





NEXT MODE WILL BE 
LOWER DUCT 
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WKB SUBROUTINE FLOW CHART 


SET UP SOME 


INITIAL 
VALUES 





= DEPTH LOOP BEGINS J = 1, NPI 

















INTEGRATE CLOSE REAL 



















IMAGINARY VERTICAL 
VERTICAL WAVE NUMBER 
WAVE NUMBER. INTEGRATION ae 
START IMAGINARY WAVE NUMBER 
WAVE NUMBER 
INTEGRATION 


CLOSE IMAGINARY 
VERTICAL 
WAVE NUMBER 
INTEGRATION 
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2 V 


COMPUTE UPPER 
PHASE ANGLE 





ARE THE 
TWO DUCTS 
ONNECTED ? 












COMPUTE PHASE 
EFFECT OF THE 
BARRIER 





RE-START INTEGRATING 
REAL VERTICAL WAVE NUMBER 





21 y 


DEPTH LOOP ENDS 


y 







BOTTOM 
REFLECTED? 


es ate, 












COMPUTE LOWER PHASE 
ANGLE FOR DECAYING 
SOLUTION 






COMPUTE BOTTOM REFLECTED 
PHASE ANGLE 






COMPUTE 
DIFFERENCE 
FUNCTION 






107 


ONE MORE DUCT 


LOWER PHASE 
= 71/4 





Memeo CCCCCOC CEC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C 
c 
NCRMO1L = A NORMAL MODE EIGENVALLE C 
AND EFGCENFUNCTION SOLVER 

e 

C 

C 


LT KIRK E. EVANS, USN 
ume Cen cCCCCOCC COG@CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


IMPLICIT REAL#8(D) 
KL KKK KKK EK SKE KE COMMON BLOCKS 35355 


CCWMON /DBLI/ DV S@;ORORB+OCKAP,DUMAX 
CCMMON Ci ers 


C 
g 
C 
e 
a 
C 
e 
C 
C 
C 
C 
C 
Cc 


ued 

— eee 

Soe 
Ze 

ate 
Om Om om 


OAZNMNZOOO 


OcnmorKrCZzOe INOCTMN 


rad 


iM 


© 


AX,OKL,DOKU 


ie Oemw ~*~ NOOFE 


ZUN ARNO SZ HUWIN 
OUw~w sO) 


o> 
on OM 


NO 
~ ROMO ZR Ze AZ NRA 


DADO QAOOMOAANIMOMS 
QRAQMAMOIMCMOAO 

i= £m. ao. ec. = 

~ AL S S SSS  S 
PNDOWMOVO—HTNHNVOTOE 
ODOLIONA SZenorawe 
OmMmAMDSZAYypvZrcecArrr 
NMA CHvCceczZzrenrmimmMm 
CN Tt 


9 
K 
1 
ZB(20) ,UB( 20) 
ant 


ee 
TIT 
THOM 
rOx< 
“Teun 


C 
CCXML KKK SKK EK KE DATA 3222999929922 >> 


DATA DP2/6.28 3185306D0/ 
DATA NPRINT/6/ 


KS KKK EEE KEK KEKE EK EK PROGRAM BEGINS 35> 


MDM COON 


WRITE (NPRINT,16) 


1 (IOUT,NFREC) 
N=DCMIN 


GO 


e 
{ 
+ 
Fei) Ore 

L €CZ22,C0CC,NP1) 


TZOYnNnNPOorYr AM 


ee a eo eee FKECUENCY LCCP EEGINS 


OOAOAGY 
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DW/DCMIN 
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N 
Zz. Fas 
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=D ™ 
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WO 
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! 
OOW 


WOU 


INT,18) FR 
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PR 
OP 
MI 


ZOOM 


MCDE LOOP BEGINS 


ene ce bees cum cee et ee ce ge ee es ce ce es es ce ee et ee es ee ee es et — = 


OUOOW 


l1»NMOD 


oe M= 
LL HALF (My, &9) 


DC 
IT 
CA 


ITERATICK LOOP EECINS 


YOU 


e5DC/DKN 


I=-DEPSIED) Uls63¢ 


ERAT (DOKN,MS) 


41898 


eer (MS-M) 


WO 


=) 


€ CKL=OKN 
GC T0 4 
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aa 
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=. 
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a 
uJ 
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© 
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| 
‘ 
1 
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Sher 
CK SKK EEE EK EEK KEKE FORMAT STATEMENTS 3D 
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NE INPUTL CTOUTsNFREQ) 


REAL*8(D) 
ce he sig a Hee oe ate ie aes st oe Se fs ak a ke oe Ae co tee ak a aie ake a oe ke ake aK a BR ob ok a Me ak a AC 3 ok oe 


TI 
at 


Vy te te Hott th 3t 


YVYOUUVIOVOVUO 


SCUND AND RUN DATA; 
CCASTANTS FCR 


iE 
QME 
He Ot Hz He aie He Hs aK oe Be SS SC See ok Ke aye he he a Neck oe he aie hc speck she ake af ak whe ie oe a He a Oe Oe ae he He ae ae soo ae ok oe OE ax 
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GQ 
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we Ver 8h VS VO eo 
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(se ees QAOW Us 
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> a ee a 
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NFTS20,RO,RE,CMIN,CB,ZMy,PLIMAX 
BZZ 


CHel2) EXTRA 
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Zh 


NPU 
fea a 
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te. 
N=OBL 
WRITE (NPUNCH, 12) 


ln 
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lL WRITE (NPUACHs,10) MODE,FQ,DK,DPVEL,OGVEL,CS,CECAY 
RETURN 





SLERCUTINE CZFLOT 
: IMPLICIT REAL*8(D) 
C De Se TK NS pe Ye Ke OE Ole oie A oe Oe i Os So oe aK ie rie ok Se ok as a eC oe oe oe ek ie a OO Oe Me Ok OE OR HK ke 
C * THIS SUBRCUTINE PLCTS THE SCLAD VELCCITY 
c * PROFILE CN THE PRINTER PLOT USING 
C # TRE N.P.S. ROUTINE UTPLCT. 
‘9 He Dis He ge a Ne os A OC Ee RK OK CK Oe Oe oe He he Me ie He ae eK oie oe Re oe aK aXe He ee Yh aE A ap Ok ae Oe mR tH 2K OK 
c 
ee <<<<<<<<<<<< KKK KEKE KKK EEK EEK KKK KKK COMMON BLECKS >>>>> 
CCMMON /CBL1/ DVZ( £01) ,0H,0H2,0CBSCE »DRORB yCKAF2,DUNAX 
COMMON /DELE2Z/ DW2,0CB,OCHIN,»OW,DC1L,ODCLSC 
COMMGN /OBLE3/ DUZ(5CL} ,OZZ(501) 
CCMMGN /PARAM/ N,NP1,IEX 
CCNMON /FEAD/ HEC(20) 
COMMON /WORK/ WRK1(20) 
COMMON /SINGL/Z ZMyCByCMIN,CMAX, PLTMAX,UM,FOQ 
CCMMON /SPLCT6/ KANGE(4) 
CCMMGN /DSCUND/ DCC(5C1) 
COMMGN /BCTTO/ DECAY, 26(20),UB( 20) 
CK KKK KE KKK 6K EEK EE EEK KK KEE KK KK KKK ERCGRAM EEGINS >>>>>> 
C 
C *%x WRITE THE HEADINGS **% 
WRITE (6,5) HED 
WRITE (636) 
: WRITE (6,7) 
IF (PLTMAX.LT.«ZM) PLTMAX=ZM 
; FE=(ZM=FLTMAX3=0.5C-1 
C ------ ~+—----- ~--+---~+- DEPTH INCREMENT LCCP BEGINS 
C 
HE 1. 1=142¢ 
WRKI(T)=CB 
ZE(1)=HBXFLOAT(I-1) 
; 1 CCNTINUE 
C -------- ~----- ar---3-- CEPTH LCCF ENDS 
C *e% COMPUTE THE RANGE FACTCRS **% 
RANGE{1L)=AMAXE(CB,CMAX) 
RANGE(2)=CMIN 
, RANGE(3)=2ZM 
C ee WILL BOTTOM PRCFILE BE DRAWN? #** 
F IF (ZM~PLTMAX) 39272 
: | 
2 MCDCUR=0 
RANGE(4)=0.0 
A GC TO 4 
2 MCCCUR=1 
: RANGE(4 )=ZM-PLTMAX 
C *#t% DRAW THE BOTTCM C(Z) PROFILE *** 
CALL UTPLOT (WRK1,Z8,20,RANGE, 1yMCOCUR) 
é MCCCUR=3 
C xe DRAW TEE OCEAN C(Z) PROFILE *** 
C 
; 4 CALL UTPLCT (CCC,CZZ,N\P1,RANGE,2,MCOCUR) 


RETURN 
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